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ABSTRACT 


Investigation  of  two-dimensional  transonic  flows  is 
extended  to  axi- symmetric  problems.  This  is  of  considerable 
practical  interest,  for  example,  with  regard  to  missiles  or 
aircraft  engines  which  approximate  much  more  closely  bodies  of 
revolution  than  two-dimensional  bodies.  The  main  concern  with 
axi -symmetric  flows  lies  not  only  with  the  complexity  of  the 
governing  nonlinear  partial  differential  equation  which  is 
mixed  of  elliptic-hyperbolic  type  but  also  with  the  lack  of  a 
general  method  for  accurately  solving  this  type  of  equation. 
We  solve  the  nonlinear  transonic  equation  using  separation 
of  variables  technique,  which  yields  two  nonlinear  ordinary 
differential  equations.  The  x-dependence  can  be  integrated 
numerically,  and  the  solution  for  the  r-dependence  can  be 
obtained  using  the  expansion  method  originated  by  Van  Dyke. 
This  works  well  with  only  three  terms  in  the  expansion.  The 
sonic  solution  of  these  equations  is  obtained  analytically 
since  both  equations  are  simple  enough  to  be  integrated  for 
this  case  (H„=1.0).  The  small  parameter  (1-Ma,^)  plays  an 
important  role  in  specifying  the  shape  of  the  boundary 
surfaces  for  external  axi-symmetric  steady  flow  of  interest 
for  design.  A  Navier- Stokes  solver  was  used  to  compute  the 
inviscid  flow  to  confirm  our  results  in  the  region  over  the 
surface  where  the  small  perturbations  apply. 
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I, 


INTRODUCTION 


Transonic  flow  is  a  classical  problem  in  gas  dynamics  and 
yet  not  many  exact  solutions  are  available.  The  main 
difficulty  lies,  of  course,  in  the  nonlinearity  of  the 
equations  and  in  the  boundary  conditions. [Ref.  1] 

It  is  of  interest  to  extend  two-dimensional  flow  solutions 
to  axi- symmetric  problems,  which  are  more  practical  for 
designing  bodies  of  revolution  such  as  missiles  or  aircraft 
engines.  The  goal  is  to  design  a  shock- free  body  or  surface  so 
we  can  minimize  the  penalty  of  the  wave  drag  that  is  caused  by 
shock  waves.  The  small  perturbation  technique  may  be 
implemented  for  simplifying  the  governing  equations  and  the 
mathematical  procedures  for  solving  the  resulting  equations  in 
nonlinear  flows.  This  research  shows  a  solution  to  nonlinear 
two-dimensional,  external  flow,  over  an  afterbody  axi- 
symmetric  surface.  Beginning  with  the  transonic  equation  for 
an  axi -symmetric  body  in  Chapter  II  and  applying  separation  of 
variables  gives  two  nonlinear  ordinary  differential  equations 
which  lead  into  two  exact  solutions.  One  of  the  nonlinear 
ordinary  differential  equation  can  be  integrated  numerically, 
and  other  one  can  be  solved  using  the  outer  expansion  method 
by  Van  Dyke  [Ref.  2]  .  Chapter  III  shows  the  derivation  of 
pressure  coefficient  and  local  Mach  nxnnber  for  different 
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surfaces.  Supersonic  boundary  surfaces  are  obtained 
implicitly  in  Chapter  IV  for  W.=1.05,l.l,  and  1.2  and  they  are 

examined  by  a  3D  Navier-Stokes  solver,  irun  as  an  Euler  solver 
for  this  axi- symmetric  problem  in  Computational  Fluid  Dynamics 
(CFD) .  In  addition,  the  sonic  surface  is  derived  analytically 
in  Chapter  V.  This  sonic  boundary  surface  is  also  examined 
using  3D  Navier-Stokes  solver  in  CFD  to  verify  the  results. 
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ZI.  TRANSONIC  EQUATION  FOR  AXI- SYMMETRIC  BODY 

The  small  perturbation,  non-linear,  axi- symmetric, 
transonic  equation  is  written  as  [Ref.  3] 

(1-Mi)<p^+ (rq>^)  =<p^(p^  (1) 

Where  the  non-dimensional  axial,  x,  and  radial,  r,  velocity 
potentials  have  been  defined  as  [Ref.  1] 

<p^  =  Ad  (y+D  ^  (Y+D  <2) 


where  is  free  stream  Mach  number  and  y  is  ratio  of  heat 
capacities . 

An  exact  solution  to  Equation  (1) ,  the  transonic  equation, 
can  be  found  from 

<p(x,r)  =<p®(x,r)  +  {1-Ad)x  (3) 


Which  patterned  after  the  two  dimensional  case  given  by 
Biblarz  [Ref.  l]  .  Here  (p®  (x,r)  is  a  separable  solution  for 
the  velocity  potential. 


3 


The  transonic  equation  in  axi- symmetric  form,  Equation  (1) , 
may  be  separated  by  j  atting  the  potential  function  <p  (x,r) 
equal  to 

<p  (x,  r)  =  5  (x)  (  (r)  +  (l-Afi)x  <<> 


Substituting  the  above  function  into  the  transonic  equation 
results  in  two  ordinary,  second  order,  non-linear  differential 
equations 


dx  dx^ 


-X  I  =0 


(5) 


dz^ 


IK 

r  dr 


-  X  C"  =  0 


(6) 


where  X  is  the  separation  constant. 

The  solution  to  the  first  O.D.E.  (5)  is  obtained  by  multiplying 
both  sides  by  d^/dx, 


kIk 

^  dx^ 


-K 

dx 


(X5)  =o 


(7) 


or 


_d 

dx 


1 

3 


[dxj  2  ^ 


(8) 
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Thus 


dx 


Rearranging 


dx  = 


di 


2 


a 


1 


and  integrating 


X-Xq  = 


(9) 


(10) 


(11) 


where  a  and  are  integration  constants. 

The  solution  to  the  second  non-linear  O.D.E. (6) ,  is  obtained 
by  an  outer  expansion  method.  (Ref.  2] 

^(r)  =— i-  +  (l-Mi)fi  (r)  +  +  .  .  .  (12) 

Xr^ 

Where  represents  a  small  parameter  and  the  first  term 

is  the  purely  sonic  solution. 
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By  taking  first  and  second  derivatives  and  substituting 
them  into  the  original  2nd  order,  non-linear,  O.D.E. (6) ,  the 
solution  becomes 


C(r) 


1 

kr^ 


4  +  |l-Ari|aAr  ^^*2) 


28+8/5 


(By 


where  a  is  constant. 


Equations  11  and  13  are  solutions  to  the  differential 
equations  5  and  6.  Boundary  conditions  are  needed  to 
determine  the  constants  a ,  ,  a,  and  k  . 

The  constants  ,  a,  and  X  can  be  shown  to  be  related  by  the 
expression  [Ref.l] 


^  C’l'l  a 


1. 08x10'' 


(14) 


and 


a  =  + 


(15) 


The  positive  sign  above  is  for  M_  *  l.O  and  the  negative  sign 
is  for  Mj,  <  1,0. 
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Upon  inserting  Eqn. (15)  into  Eqn. (11) ,  we  have 


ci^ 


(16) 


Factoring  out  ,  we  obtain 


x-x^ 


_ _ 


(17) 


Introducing  new  variables 


I  =  I 


(18) 


^(v^+2) 


(19) 


Equations  (17)  and  (13)  become 


(20) 
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+2)  ^  ^  ^  ^  ) 


[  '  '  50.63 

Finally,  defining  two  new  variables 

^  ^  '  -  (4)" 


(a  X) 


Equations  20  and  21  become 


‘=  f 


where  Xq=0. 


2  2 

?  =  4  +  h-M^I  +iU^lL  fZCv/ff.z) 

^2]  I  '  50,63 


Equation  (24)  will  be  numerically  integrated  and  plotted  in 
Fig.  (1)  as  I  (^)  vs.  x  ,  and  Equation  (25)  will  be  evaluated 
and  plotted  in  Fig. (2)  as  ?  (f)  vs.  f  for  M.=  1.05,  l.l, 

1.2.  These  results  will  be  used  later  to  obtain  transonic 
surfaces  with  their  corresponding  C  and  local  Mach  number 
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profiles.  In  Fig.  (2)  ,  Equation  {25)  is  plotted  up  to  its 
minimum  value  with  increasing  f.  Thereafter,  a  constant 
value  is  patched  in  to  the  right  of  the  +  mark.  This  is  done 
because  of  the  need  to  keep  C  constant  reflecting  the 
necessary  behavior  of  the  function  for  range  f  [Ref .1] . 
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III. 


PRESSX7RE  COEFFICIENT  AND  LOCAL  MACH  NUMBER 


The  pressure  coefficient  for  a 
perturbations  is  given  by  [Ref.  4] 


flow  with  small 


(1-Ad) 


(26) 


The  linearized  pressure  approximation  for  axi- symmetric  flow 
is 


-2 

Urn 


(27) 


Recall  the  axial  velocity  potential.  Equation  (2)  ,  Chapter  II, 


<Px  =  ^  (Y+1) 


(28) 


thus 


u 


9. 


Ad  (y+1) 


(29) 
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siibstitute  Equation  (29)  into  Equation  (27)  ,  yields 


-2  <Px 

Ml  (y+1) 


(30) 


We  need  to  take  the  derivative  of  the  potential  function  with 
respect  to  x  in  equation  (4) . 


<Px  = 


c 


dx 


(31) 


Rewriting  equation  (9)  with  the  constant  a  inserted  becomes 


dx 


3X 

2 


e  +  q  !  (1-M1)  |l-7574l3 


(32) 


Recall  Equation  (10) 


I  =  I 


3A  ' 

2CJ 


(33) 


Inserting  Eqn. (33)  into  Eqn. (32)  and  factoring  out,  yields 


_  A 

~dJc  - 


[  V  +  \a-Ml)  |i-"574  ] 


(34) 
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Recall  Equation  (23)  and  rearrange 


(a 


(35) 


Rewriting  Equation  (31)  the  potential  function  as 


C  (1-^)  ^  +  (i-Ar«) 


Thus 


<Px 


J. 

30.4142^^3 

j^FTiais 


(37) 


Recall  Equation  (14)  Chapter  II 


^  3  1.2426 


1.08x10'2 


(38) 


then  finally  Equation  (37)  becomes 


<p^  =  0.2208  C  [  +  |(l-Mi) 


(39) 


Rewriting  the  coefficient  of  pressure 


C  = 


-2 


^  td{y+i) 


0,2208j[r+  |(l-wf)  p  ’'5''^]3  +  (l-wi)  (40) 


The  previous  expression  for  Cp  will  be  used  further  in 

determining  the  local  Mach  number. 

The  perturbation  velocities  are  assumed  to  be  vary  small 
compared  to  the  undisturbed  uniform  velocity  U^.  [Ref.  5] 

Hence  at  a  point  near  the  body,  the  velocity  vector  V  for 
2-D  flow  is  given  by 

V  =  i{u^  +  +  j  Uj. 


We  can  obtain  the  local  Mach  nximber  M  in  terms  of  the 
perturbation  velocities  and  the  local  speed  of  sound  c  as 

(42) 


Neglecting  the  higher  order  ratios  of  the  free  stream  to 
perturbation  velocities,  since  they  are  considerably  smaller 
than  unity,  we  have 
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The  ratio  of  the  local  speed  of  sound  to  the  undisturbed 
uniform  speed  of  sound  becomes 


(44) 


or 


£!  = 

cl  '  3^0/r 


(45) 


Where  is  the  stagnation  temperature. 

We  can  write  temperature  ratios  as  functions  of  Mach  numbers 
as 


1  + 

1  + 

2 


(46) 


Substituting  Equation  (45)  into  Equation  (46)  and  rearrange, 
we  have 


m2  =  M« 


I  A  2 


(47) 
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from  Equation  (27) 


-2  u. 

C7, 


Rearranging  and  solving  for  local  Mach  number 

j,2  =  ^  (1-gp) 

1  *  Mi  Cp 

This  equation  will  be  used  to  show  the  local  Mach 
each  surface  at  M.  =  1.05,  1.1  and  1.2  [Ref.  6] . 


(48) 


(49) 


number  for 
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IV. 


SUPERSONIC  BOUNDARY  SURFACES 


The  boundary  conditions  which  need  to  be  applied  require 
that  the  gradient  of  <p  vanish  far  ahead  of  the  body  and  that 
the  flow  be  tangential  to  the  surface  [  Ref.  7] .  In  terms  of 
perturbation  velocities  the  boundary  conditions  become 


(50) 


Recall  the  modified  velocity  perturbation  potential  Eqn. (2) 


<p,  =  fd  (y  +  1)-^ 


(51) 


then ,  we  have 

Ml  (Y+1) 


(52) 


by  substituting  Equation  (52)  into  Equation  (50) ,  we  obtain 


1^]  = 

\dxl  Ml  (Y  +  1) 


(53) 
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By  differentiating  Equation  (4)  with  respect  to  r  gives 


<Pr 


=  I  K 
^  dr 


(54) 


Recalling  Eqn. (18)  and  substituting  Eqn. (21)  into  Eqn. (23) , 
and  taking  the  derivative  with  respect  to  f ,  the  above 
equation  becomes 

=  0.0848  I  (55) 


Substituting  Equation  (55)  and  (13)  into  equation  (53),  we 
have 


0.0848  ^  ^8 
Alf(Y+l) 

+  0,1512  (1 


+  2.8284  ll-A^l 
-Mf)2  f®-®” 


^1.8284 


(56) 


It  can  be  noticed  here  that  Z  and  Z  are  given  implicitly  as 
function  of  x  and  f  respectively. 

From  Eqns. (19)  and  (22),  we  can  rearrange  Equation  (56)  as 


dSc 


—  I  ^  *  2.8284  \1-M^\  f  1-8264 

wf(Y+l)  Lf" 


+  0.1512  (1-Mi)2  f6.657 


(57) 
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Further  arrangment  of  the  previous  equation  yields 


_ df _ 

-^-2.8284|l-iy:f|  _  0 . 1512  (1-J^)  ^  _f6.657 

\5oj 

I  dx 

Mi(Y+i) 

Now,  we  can  compute  the  exact  boundary  surfaces  out  of  the 
equation  above.  Starting  with  the  left  hand  side  of  the 
equation 


Z(i')  = 


8-2 . 8284  [I-Wf  |i'^-®2®^-0.1512  (l-i^f 


(59) 


By  defining  equation  (59)  we  can  plot  Z  (f)  vs.  f  in  Fig. 

(3)  for  M.  =  1.05,  1.1,  and  1.2. 

The  minimum  values  of  different  Mach  numbers  can  be 

expressed  as  {MJ  which  has  been  found  to  be  [Ref.  1] 


1.2074 

■  1.^2  |o. 20*71 


(60) 
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Fig.  (4)  represents  vs.  W.  and  shows  a  symmetry  close  to 
/f.=1.0  where  So  for  our  purpose  to  represent  the 

condition  sufficiently  close  to  W.  =  l.O,  we  will  choose  a 
finite  value  of  {perhaps  as  4.0). 

Defining  the  right  hand  side  of  Equation  (58)  as 

K{x)  =  -3.26x10"^  r 


ta)cing  the  derivative  of  Equation  (24)  and  substitute  it  into 
Equation  (61) ,  yields 


K(Jt)  = 


-3 .26x10'^ 
Mi  (y+1) 


Ul 


(62) 


integrating  by  parts,  we  obtain 


K{x)  = 


-2.45x10'^ 

(Y 


+  1)A1^  / 


1716 


(€3) 


the  equation  above  will  allow  us  to  plot  K{ic)vs.^  in 
Fig. (5)  for  ■  1.05,  1.1,  and  1.2.  Where  K (x) <  0 

for  M,*  1.0  and  K  (x)  >  0  for  <1.0. 
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Therefore,  integrating  both  sides  of  equation  (58)  becomes 


_ df _ 

—  ~2 .828411 -Mi 1512  (1 -Mi) 


-3.26x10'^ 

Ml  (y+1) 


I  dX 


(64) 


Equation (64)  is  numerically  integrated  and  which  then  allows 
us  to  determine  the  boundary  surfaces  in  dimens ionalize  and 
non- dimensional ize  (normalize)  form  as  shown  in  Figs.  6  and  7 
for  =  1.05,  1.1,  and  1.2.  Now,  we  can  determine  Cp  from 

Eqn.  (40)  and  local  Mach  number  from  Eqn.  (49)  for  the  three 
surfaces  obtained  at  different  Mach  numbers. 

In  Figs.  8-13,  the  supersonic  boundary  surfaces  for 
t^=l . 05 , 1 . 10 ,  and  1.20  are  plotted  with  their  pressure 
coefficients  and  local  Mach  profiles  in  dimensionalized  and 
non-dimensionalized  (normalized)  form.  These  results  will  be 
verified  later  using  a  3-D  Navier-Stokes  solver  in 
Computational  Fluid  Dynamics  (CFD) . 
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V. 


SONIC  BOUNDARY  SURFACE 


The  sonic  flow  (Af^=1.0)  can  be  derived  in  explicit  form 

from  Equations  (24)  and  (25) .  The  resulting  equations  for 
sonic  flow  are 


Jt  = 


(65) 


and 


(66) 


Rearranging  Equation  (65) 


Jf  = 


2 

3  dz 


(67) 


since  Equation  (67)  is  integrable,  we  obtain 

X  =  2  V 

or 


(68) 


(69) 
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We  can  write  Equations  (11)  and  (13)  in  sonic  form  as 


X 


and 


C  = 


4 


(70) 


(71) 


integrating  Equation  (70) ,  we  have 


5(x) 


(72) 


From  Equation  (4)  ,  we  can  write  the  sonic  perturbation 
potential  as 

9  (x,r)  =  ^(x)C  U)  (73) 

Substituting  Equations  (71)  and  (72)  into  Equation  (73)  yields 

<p(x,r)  = 

Equation  (57)  can  be  written  for  W',=l .  0  as 

^  -  0.2608  I 
cK  (y  +  1) 


(74) 


(75) 
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Rearranging  Equation  (75) ,  we  have 


df  =-0.1087  I  dx 


Therefore,  integrating  both  sides  of  Equation  (76) ,  we  obtain 


dz  =  -  0.1087  I  C  dSt 


Inserting  Equation  (69)  into  Equation  (77)  yields 


=  -  0.0040ic^ 


where  i'o=4.0  has  been  chosen  as  an  asymptotic  value  to 

represent  the  condition  sufficiently  well  at  W,=1.0.  Equation 

(78)  will  allow  us  to  determine  sonic  boundary  surface. 

The  pressure  coefficient  can  be  written  from  Eqn. (40)  as 


-  _  -  0 .4416  y  y  3 
(Y+1) 


in  case  of  W.=l . 0 . 
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The  local  Mach  number  can  be  written  from  Equation  (49)  for 
the  sonic  surface  as 


W2  =. 


-  Q 


1  -*■ 


111 

2 


(80) 


The  sonic  surface  is  represented  dimensionally  in  Fig.  (14) 
with  the  pressure  coefficient  and  local  Mach  number  profiles, 
and  normalized  in  Fig.  (15)  with  ^^=4.0  for  (Mjo=1.0)  .  The 

sonic  boundary  surfaces,  in  both  forms,  are  plotted  with  the 
supersonic  boundary  surfaces  in  Figs.  16  and  17. 
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VI. 


COMPUTATIOKAL  FLUID  DYNAMICS  (CPD) 


The  history  of  Computational  Fluid  Dynamics  (CFD)  is 
closely  tied  to  the  rapid  advances  in  digital  computer  which 
has  a  great  impact  on  problems  of  design  in  modern  engineering 
practice.  These  problems  can  now  be  solved  at  very  little  cost 
in  a  few  seconds  of  computer  time  which  would  have  taken  years 
to  work  out  with  the  computational  methods  and  computers 
available  twenty  years  ago.  [Ref.  8] 

The  CFD  will  be  used  to  compute  the  axi- symmetric  flow 
over  the  boundary  surfaces  obtained  by  the  small  perturbation 
method . 

The  computer  program  GRAPE  [Ref .9] ,  an  acronym  derived 
from  Grids  about  Airfoils  using  Poisson's  Equation,  has  been 
written  by  R.  Sorenson  at  Ames  Research  Center  to  generate 
two-dimensional  finite  difference  grids  about  airfoils  and 
other  shapes  by  the  use  of  the  Poisson  differential  equation. 
Outer  and  inner  boundaries  are  specified  for  the  C-type  grid, 
where  the  surface  of  the  body  is  treated  as  the  inner 
boundary.  Important  characteristics  in  grid  generation  are 
the  ability  to  specify  the  spacing  between  mesh  points  at  the 
boundary,  in  the  direction  normal  to  the  boundary,  and  the 
control  of  the  angles  with  which  mesh  lines  intersect  the 
boundaries  which  is  known  as  orthogonality. 
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The  grid  size  for  the  sonic  surface  (M„=1.0)  is  107  x  60  and 
for  the  supersonic  surfaces  (M„=l.l,  1.2)  is  115  x  60  as 
shown  in  Fig.  (18)  .  Since  the  surfaces  are  symmetrical,  half  of 
the  grid  or  eventually  the  lower  surface  was  rotated  for  11 
planes  plus  two  more  for  symmetry  to  generate  an  axi- 
symmetric  after  body  surface. 

The  OVERFLOW  program  [Ref  .10]  developed  by  Ames  Research 
Center  which  uses  3-D  Navier-Stolces  and  Euler  solver  for 
viscous/  inviscid  flow  was  applied.  The  results  are  shown  in 
Figs.  19-21  for  three  different  boundary  surfaces  at  M*  =  1-0, 
1.1,  and  1.2  with  1500  iterations.  It  was  found  that  the 
density  residuals  decreased  by  more  than  one  order  of 
magnitude  over  1500  iterations.  After  3000  iterations,  the 
solution  converged  by  two  orders  of  magnitude. 

In  case  of  the  sonic  surface  Fig.  (19)  with  M,„=1.0,  the 
Mach  lines  contour  shows  the  shock  is  forming  downstream  at  a 
local  Mach  M=1.65  at  which  the  flow  becomes  subsonic 
downstream  of  the  shock.  This  result  complies  with  the  small 
perturbation  transonic  solution  obtained  earlier  in  Fig. (15) . 
In  other  words,  the  flow  is  shock  free  over  the  sonic  surface 
in  for  r^g,»1.0. 

For  the  supersonic  boundary  surfaces  M5o=1.10,  and  1.20,  we 
can  see  that  the  flow  starts  at  M=1.10  and  M=1.20  respectively 
and  follows  the  afterbody  surface  until  it  shocks.  The  shock 
is  formed  at  local  Mach  M=1.75  and  1.85  respectively  as  shown 
in  Figs.  20  and  21. 
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On  the  other  hand,  these  results  agree  with  the  small 
perturbation  solution  over  the  transonic  range  as  plotted  in 
Figs.  11  and  13  where  the  local  Mach  is  represented  by  the 
dash  line  and  it  is  increasing  towards  the  afterbody  surface 
until  it  shocks.  This  indicates  that  the  flow  is  shockless 
over  these  supersonic  boundary  surfaces  in  the  transonic 
regime. 

The  velocity  vectors  are  plotted  in  Fig.  (22)  for  Meo=1.0 
which  shows  that  the  flow  is  following  the  afterbody  surface 
until  the  shock,  and  then  flow  separation  takes  place  due  to 
the  steep  pressure  rise  and  the  introduction  of  numerical 
viscosity  into  the  flowfield  at  this  location. 

These  boundary  surfaces  obtained  are  of  interest  for 
design  of  body  of  revolution  such  as  missiles  afterbody  or 
aircraft  engines,  in  which  a  patching  technique  may  be  applied 
to  these  surfaces  to  get  the  best  shockless  surface  in  the 
transonic  range. 
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VII. 


CONCLUSION  AND  RECOMMENDATION 


We  have  shown  in  this  research  a  solution  to  the  non¬ 
linear  transonic  small  perturbation  equation  by  implementing 
the  separation  of  variables  technique.  The  x- dependence  was 
integrated  niimerically  and  the  r- dependence  was  solved  by  an 
outer  expansion  method.  It  was  found  that  the  parameter  (l- 
has  strong  effect  in  specifying  the  shape  of  the  boundary 
surfaces  of  interest  for  design.  A  3-D  viscous/  inviscid 
Navier- Stokes/  Euler  solver  in  Computational  Fluid  Dynamics 
(CFD)  confirmed  our  results  obtained  from  the  small 
perturbation  technique  for  the  boundary  surfaces  at 
1.10,  and  1.20  .  In  other  words,  we  achieved  our  goal  of 
designing  shockless  surface  in  the  transonic  range. 

This  research  may  be  extended  using  a  patching  technique 
to  obtain  the  best  shockless  surface  in  the  transonic  flows. 
Patching  criteria  will  have  to  be  developed  based  on  the 
mathematics  and  the  flow  constraints.  Also,  connection  between 
small  perturbation  solution  and  CFD  might  eventually  be  done 
internally  in  the  computer.  The  small  perturbation  is  the  most 
efficient  method  to  define  or  to  design  axi- symmetric 
afterbody  transonic  surfaces,  however  CFD  can  be  used  to 
predict  the  performance  of  these  aerodynamic  surfaces. 
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Zeta(?)  vs.  r  For  M.=  1.05,1.10,1.20 


Figur*  2.  Munarical  solution  of  Egn.  25. 
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Pigura  4.  Nunarical  solution  of  Eqn.  60 


K(!0  vs.  X  For  M„=  1.05,1.10,1.20 
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Figure  5.  Nuaerical  solution  of  Eqn.  63. 
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Figur*  6.  Numerical  solution  of  Eqn.  64. 
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Figure  7.  Nuacrical  solution  of  Eqn.  €4  (nornalisod) . 


Cp  &  local  Mach  FOR  M.=  1.05 
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Figur*  8.  cp  and  local  Mach  for 


Cp  &  local  Mach  FOR 
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Flgur*  10.  Cp  and  local  Mach  for  M,.  s  1.10 


Cp  &  local  Mach  FOR  Ma,=  1.1  ,  ro  =  1.669 


10  (nornalixad) 


Cp  &  local  Mach  FOR  M.=  1.20 


Figur*  12.  Cp  and  local  Maeb  for 


Cp  &  local  Mach  FOR  1.20  ,%=  1.431 
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Figurs  13.  cp  and  local  Mach  for  s  1.20  (normalized) 


Cp  &  local  Mach  FOR 
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Figur*  14.  Cp  and  local  Hach  for 


Cp  &  local  Mach  FOR  M^=  1.0  ,  ?o  =  4.0 
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Figur*  15.  Cp  and  local  Mach  for  M,  =  i.O  (normaliiad) 


BOUNDARY  SURFACES 


6 


Flgura  17.  Boundary  surfacas  for  M.  =  i.o,1.05,l.l»1.2  (normallzad) . 
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Figure  18.  Afterbody  axl-syimnetrlc  grid. 


FIGURE  n.  HACH  CONTOURS  FOR 


FIGURE  aa.  HACH  CONTOURS  FOR  N  =1.10 


I 


I 
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FIGURE  51.  flACH  CONTOURS  FOR  M  =1.50. 
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Figure  22.  Velocity  vectors  for 


APPENDIX  B 
COMPUTER  PROGRAMS 


*  * 

*  This  program  is  designed  to  calculate  KSI (X) ,  * 

*  K(X),  and  X  using  numerical  integration  based  * 

*  on  trapezoidal  rule  to  solve  Eqn. (24)  and  * 

*  plotting  the  output  as  shovm  in  Figs.l  &  3  .  * 

*  ♦ 

PROGRAM  KSI 

REAL  M(3) ,X(0:401,3) ,A,A1,B(401,3) 

+  ,KSI,P,FUNC,H,DD,Q,L,K(401,3) ,N,N1 

INTEGER  YY 

OPEN (UNIT=9 , FILE* ' KSI ' , STATUS* ' UNKNOWN' ) 

PRINT  *, 'ENTER  LOWER  fc  UPPER  BOUND,  #  OF  INTERVALS, 

+#0F  DATA  SET' 

READ  *,  A,  Al,  N,  N1 

100  DO  10  1*1,3 

PRINT  * , ' ENTER  MACH  NO .  ' 

READ  ♦,  M(I) 

C  M(I)  *  0. 9+1*0. 1 

IF  (M(I) .LT.1.0)  THEN 
P*-l. 

ELSE 

P=l. 

ENDIF 

DO  20  J=1,N1 
B(J,I)  *  J*A1/N1 
H  *  (B(J,I)-A)/N 
AREA  «  0. 

K(J,I)-(-0.0102/(M(I)**2))*((ABS( (B(J,I)**2) - 
+( (ABS (1-M(I)**2)) **1.7574)))** (2. /3.)+ 

+ (ABS(1-M (I) **2) ) **1.1716) 

DO  30  L-1  ,  N 
KSI  -  A  +  (L-0.5) *H 

DD  *  (KSI**2+P* (ABS (l-M(I) **2) ) **1.7574) 

IF  (DD.GT.O.)  THEN 
FUNC  *  l./(DD**{l./3.)) 

Q  -  1. 

ELSE 
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FUNC  »  1./((-1*DD)**(1./3.)) 

Q  ■■  - 1 . 

ENDIF 

AREA  =  AREA  +  H*FUNC*Q 
X(J,I>  -  AREA 

30  CONTINUE 

20  CONTINUE 

10  CONTINUE 

DO  40  J-1,N1 

PRINT  50,  (B(J,I) ,X(J,I) ,K(J,I) ,1-1,3) 

50  F0RMAT{1X,3 {1F6.3,2F10.5,1X) } 

40  CONTINUE 

200  PRINT*, 'TRY  AGAIN  ?  ENTER  1  FOR  YES, 2  FOR  DATA 

+FILE, OTHERS  FOR  NO' 

READ  *,  YY 
IF  (YY.EQ.l)  THEN 
GO  TO  100 
ELSE 

IF  (YY.EQ.2)  THEN 
DO  110  J=1,N1 

WRITE  (9,50)  (B(J,I) ,X(J,I) ,K(J,I) ,1=1,3) 

110  CONTINUE 

GO  TO  200 
ENDIF 
ENDIF 
END 
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**************************************************************** 
*  * 

*  This  progreun  is  written  to  calculate  ZETA(r)  and  ♦ 

*  r  using  Eqn. (25)  and  plotting  the  output  in  Fig. 3  .  * 

*  * 
**************************************************************** 


program  zeta 

real  M{3) , zeta (0 : 14 , 3) , r (14) 
open(unit*9,  file='  zz' ,  status* 'unJcnown' ) 

do  10  1=1,3 

c  M(I) =0. 7+1*0. 1 

print  * ,  ' enter  mach  no . ' 
read  * ,  M ( I ) 

do  20  J=l,14 
zeta (0 , I) =1000 
print  * , ' enter  r  values ' 
read  *,  r(J) 
c  r(J)=0.1*J 

zeta(J,I)=(4/(r(J)**2) )+  (abs(l- (M(I)**2) )* 

+(r(J) **2.8284) )+((!- (M(I)**2) ) **2) * (r ( J) **7 . 657) /50 . 63 
if  (2eta(J,I) .GE.zeta(J-l,I) )  zeta (J, I) =zeta (J-l, I) 

20  continue 

10  continue 

print  30,  M 

30  format (llx, 5F10 .4) 

do  40  J=l,14 

c  print  50,  r (J) , (zeta (J, I) , 1=1, 4) 

write  (9,50)  (r (J) , zeta ( J, I) , 1=1, 3) 

50  format  (4x, 3 (2F10 .4 , 2x) ) 

40  continue 

end 
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*  * 

*  This  program  is  written  to  calculate  Z(r)  and  r  * 

*  using  Eqn. (59)  and  showing  the  output  in  Fig. 4  .  * 

*  * 


PROGRAM  ZR 

REAL  M(4) ,R(1000) ,2(0:1000,4) 

0PEN(UNIT»99,FILE='ZR1' , STATUS =' UNKNOWN' ) 

DO  10  I»l,4 

C  M(I)  =  0.7  +  1*0.1 

PRINT  *,  'ENTER  MACH  NO.  ' 

READ  *,  M(I) 

DO  20  J=l,1000 
Z(0,I)  =  0.0 
R{J)  *=  0.001*J 

Z{J,I)  =  (R{J)**3)/{8 

(2.8284* (R(J) **4.8284) *ABS (1- (M(I) **2) ) ) - 

+(0.1512* (R(J) **9.657) * ( (l-M(I) **2) **2) ) ) 

IF  (Z(J,I) .LE.0.0)  Z(J,I)=Z(J-1,I) 

20  CONTIJJUE 
10  CONTINUE 


DO  40  J=l,1000 

C  PRINT  50,  R(J)  ,  (Z(J,I) ,1=1,4) 

WRITE  (99,50)  R(J) ,  (Z(J,I),  1=1,4) 

50  FORMAT(4X,5F10.6) 

40  CONTINUE 

END 
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*  * 

*  This  program  is  written  to  determine  Rmin  for  * 

*  transonic  Mach  numbers  using  Eqn. (60)  and  the  ♦ 

*  output  is  shown  in  Fig. 5  .  * 

*  * 


PROGRAM  RMIN 
REAL  M,R 

OPEN (UNIT=88,FILE=' RMIN'  , STATUS  =' un)cnown'  ) 

DO  100  M  =  0.8,1,2,0.0001 
IF(M.LT.l.O)  THEN 

R=1.2074/ ( (1-M**2) **0.2071) 

ELSE 

IF(M.GT.l.O)  THEN 

R=1.2074/ ( (M**2-l) **0.2071) 

ELSE 

R  »  7.0 
ENDIF 
ENDIF 

C  PRINT  80,  M,R 

WRITE  (88,80)  M,R 
80  FORMAT(1X,F7.5,7X,1F10.6) 

100  CONTINUE 
END 
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*  * 

*  This  program  is  designed  to  niomerically  integrate  * 

*  Z(r)  ,Eqn.(64),  based  on  trapezoidal  rule  .  * 

*  * 


PROGRAM  IZR 

REAL  M{3) ,R{401,3) , A, A1 , IZ (0 : 401, 3 ) , H, FUNC, N,N1 , Ro 
INTEGER  YY 

OPEN(UNIT=55,FILE='IZR' , STATUS® 'UNKNOWN' } 

100  DO  10  1=1,3 

C  M(I)  »  0.7  +  1*0.1 

PRINT  *,  'ENTER  MACH  NO.  ' 

READ  *,  M(I) 

PRINT  *,  'DETERMINE  THE  VALUES  OF  MACH  #',M(I) 

PRINT  *,  'ENTER  Rmin  ,A' 

READ  ♦,  A 

PRINT  *,  'ENTER  THE  UPPER  LIMIT  OF  R  ,A1' 

READ  *,  A1 

PRINT  ♦,  'ENTER  #  OF  INTERVALS  ,  N' 

READ  *,N 

PRINT  *,  'HOW  MANY  DATA  SET  DERIVED  ?  Nl' 

READ  *,  Nl 

PRINT  *,  'INPUT  Ro  FOR  THIS  MACH  #  ,  Ro' 

READ  *,  Ro 

DO  20  J=1,N1 
R(J,I)  =  J*A1/N1 
H  =  (R(J,I)-A)/N 
AREA  =  0. 

DO  30  L=1,N 
RR  =  A+(L-0.5)*H 

FUNC= (RR**3) / (8- (2.8284* (RR**4 . 8284 ) *ABS { 1 - (M(I) **2) ) ) 
+  - (0.1512* (RR**9.657) * ( (l-M(I) **2. ) **2. ) ) ) 

AREA  =  AREA  +  H*FUNC 
IZ(J,I)  =  Ro  -  AREA 

30  CONTINUE 

20  CONTINUE 
10  CONTINUE 

DO  40  J=1,N1 

C  PRINT  50,  (R(J,I) ,IZ(J,I) ,1-1,3) 
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WRITE{55,50)  {R(J,I) ,IZ(J,I) ,1=1,3) 

50  FORMAT  ( IX, 3 {F6 . 3 , F9 . 6 , 3X) ) 

40  CONTINUE 

PRINT  * , ' TRY  AGAIN  ?  ENTER  1  FOR  #YES# ,  OTHERS  FOR  #NO# ' 
READ  *,  YY 

IF  (YY.EQ.l)  GO  TO  100 
END 
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it 

*  This  program  is  using  numerical  integration  to 

*  determine  the  supersonic  boundary  surfaces  , 

*  Eqn. (64) ,  and  calculating  pressure  coefficient 

*  and  local  Mach  niimber  profiles,  using  Eqns.  (40 

*  &  49) .  The  output  is  shown  in  Figs. (6-13)  . 

* 


PROGRAM  CPP 

REAL  X(401. 9) ,Y(401,6) ,  XR (401) , YR (401) ,M, MX, XYMAT 
INTEGER  I, J,K,L, IWRITE 

OPEN (UNIT=52 , FILE= ' cpl5.dat ' , STATUS= ' UNKNOWN' ) 
OPEN(UNIT=53,FILE=' Cpll.dat' , STATUS UNKNOWN' ) 

OPEN (UNIT=54,FILE=' cpl2.dat' , STATUS =' UNKNOWN' ) 
OPEN(UNIT=77,FILE='cpp.dat' , STATUS =' UNKNOWN' ) 
OPEN(UNIT=78,FILE='ml05.dat' , STATUS =' UNKNOWN' ) 
OPEN(UNIT=79,FILE=' tirel.dat' , STATUS =' UNKNOWN' ) 
OPEN (UNIT=80,FILE=' tire2.dat' ,STATUS=' UNKNOWN' ) 

100  OPEN (UNIT=1,FILE='KSI. DAT' , STATUS= ' OLD' ) 

OPEN (UNIT=2,FILE='IZR. DAT' , STATUS= ' OLD' ) 

REWIND (UNIT=1) 

REWIND (UNIT=2 ) 

IWRITE  =78 
MACP=52 

READ  (1,*,END=50)  ( (X ( I , J) , J=1 , 9 ) , 1=1 , 401 ) 

50  READ  (2,*,END=60)  ( (Y (K, L) , L=l, 6) , K=l, 401) 

60  DO  10  J=2,6,2 

PRINT  *,'  Ro  is  rmin  ;  EPSILON  is  the  accuracy  ' 
PRINT  *, 'ENTER  MACH  NO.,  Ro  ,  EPSILON' 

READ  *,  M,RO,EPS 
WRITE  (77,*)  M,RO,EPS 
XM  =  1-M**2 
L  =  J/2*3 

DO  30  I  =  1,401 
DO  20  K  =  1,401 

IF(Y(I, J) .LT. 0.001)  THEN 
GO  TO  10 
ELSE 

XYMAT  =  ABS(Y(I,  J) -ABS(X(K,L)  )  ) 

IF (XYMAT.GT.EPS)  THEN 
GO  TO  20 
ELSE 


★ 
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XR(K)  -  X(K,L-1)/R0 
YR(I)  -  Y(I,J-1)/R0 

ZETA  -  (4./ (Y(I, J-1) **2) )  +  ABS(XM)*Y(I,J-1)**2 
+  +  XM**2* (Y(I, J-1) **7.657) /50. 63 

Cp  «  (-2/ (M**2*2.4) ) * (0.2208*ZETA* (X(K,L-2) **2 
+  (ABS(XM) ) **1.7574) **(l./3.)+XM) 

MX  »  SQRT(ABS( (M**2* (1-Cp) / (l+0.2*M**2*Cp) ) ) ) 

C  PRINT  90,  X(K,L-1) ,Y{I, J-1) ,Cp,MX 

WRITE(77,40)  X(K,L-1) ,Y(I,J-1) ,Cp,MX 
WRITE (MACP, 40)  XR(K)  ,YR(I)  , Cp  ,MX 
WRITE (IWRITE, 70)  XR(K)  ,  YR(I) 

GO  TO  30 
ENDIF 
ENDIF 

20  CONTINUE 

30  CONTINUE 

IWRITE=IWRITE+1 
MACP  «»MACP+1 
10  CONTINUE 

40  FORMAT  ( IX, 4 (F15 . 7 , 2X) ) 

70  FORMAT  {1X,2F17.9) 

90  FORMAT  { IX, 4 (F15 . 7 , 2X) ) 

PRINT  * , ' TRY  AGAIN  ?  1  FOR  YES ! , 2  FOR  NO ! ! ' 

READ  *,  lY 

IF(IY.EQ.l)  GO  TO  100 
END 


.8284 
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*  This  program  is  written  to  determine  the  sonic 

*  boundary  surface  M  =  l.O,  ro  =  4.0  using  Eqn. (78) , 

*  and  calculating  pressure  coefficient  and  local 

*  Mach  profile  using  Eqns.(79  &  80).  The  output  is 

*  shown  in  Figs.  16  &  17  .  Some  output  files  are  used 

*  in  CPD. 


★ 
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program  sonic 

real  M(0:160) ,r(0:160) ,x{0:160) ,zeta (0:160) ,ksi (0:160)  , 
+  2(0:160) ,k(0:160) ,Cp(0:160) ,  RR ( 0 : 160) , YY ( 0 : 160) , 

+  XX (0 : 160) , XXX (0 ; 160) 


c 

c 

c 

c 


c 

10 


open (unit=ll, f ile= ' sonicl.dat' , status= ' unknown' ) 
open (unit=22 , f ile= ' sonic2.dat' , status^ 'unknown' ) 
open(unit=33 , f ile=' sonic3.dat' , status=' unknown' ) 
open (uni t =66, file=' sonic4.dat' , status= ' unknown' ) 


sonicl.dat 
sonic2 . dat 
sonic3 .dat 
sonic4.dat 


=>  X  ,  r  ,  Cp  ,  M 
=>  X  ,  r  (horintal) 

=>  X  ,  r  (vertical) 

=>  x/ro  ,  r/ro  ,  Cp  ,  M 


do  10  I  =  0,  159 
x(I)  =  1*0.1 
ro  =  4.0 

r(I)  =  (abs(ro**4.  -  0.004* (x(I) **4. ))) **0.25 

RR(I)  =  r(I)/ro 

YY(I)  =  -  1.0*  RR(I) 

xx(I)  =  x(I) /ro 

xxx(I)  =  3.975  -  xx(I) 

zetad)  =  4./(r(I)**2.) 

ksi(I)  =  {x(I)**3.)/27.0 

2(1)  =  r(I)**3. 

k(I)  =  -  0.00101* (x(I) **4. ) 

Cpd)  =  -  0.1840  *  zetad)  *  (ksid)  **0.6667) 
M(I)  =  (abs( (l-Cp(I))/(l  +  0.2*Cp(I) ) ) ) **0.50 
print  *,  xxd)  ,RR(I)  ,Cpd)  ,M(I) 
continue 


do  20  I  =0,  159 

write(ll,50)  xx(I)  ,  rr(I)  ,  Cp(I)  ,  M(I) 
write(22,60)  xxx(I)  ,  RR(I)  ,  YY(I) 
write(33,70)  xxxd),YY(I) 
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write(66,5C)  x(I)  ,r(I)  ,Cp(I)  ,M(I) 

50  format {2x, If 8 . 4 , 2x, 3f 14 . 7, 2x) 

60  format (2x, lf8 . 3 , 2x, 2f 12 . 6, 2x) 

70  format (lOx, If 12 . 6 , 5x, f 12 . 6) 

20  continue 

write(34,80) (xxx(I) ,i=0,159,3) 
write (34, 80) (xxx(i) , i=156, 0, -3) 
write (34,80) (yy(i),i=0, 159 , 3 ) 
write (34,80) (rr(i),i=156,0,-3) 

80  format (5 (lx, fl2 . 6, ' , ' ) ) 

do  30  I  =158,  0  ,  -1 
write(33,70)  xxx(I),RR(I) 

30  continue 

end 
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